{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Author:__ Bram Van de Sande\n",
    "\n",
    "__Date:__ 25 APR 2018\n",
    "\n",
    "__Remaining implementation challenges:__\n",
    "    \n",
    "1. Disk-volume of region-based databases: the inverted database design can significantly reduce the size on disk (from 120Gb to 4,7Gb for the 1M regions-24K features human database) but with a huge impact on read performance (from 1 second to several minutes for a typical signature). Potential mitigation challenges are:\n",
    "    - Recoding the decompression of the inverted design using C++ and its standard template library. Impact only be a constant factor while personal investment will be substantial.\n",
    "    - Break the clean interface between database storage and AUC and recovery curve calculation. The exact order of the genes is not really required for AUC calculation so decompression times could be significantly reduced.\n",
    "2. Large memory footptint when using many cores combined with large rank_threshold increase the probability of generating MemoryError's. This can be mitigated by reducing the number of cores used (current strategy). The burden on memory can be reliefed by reimplementing the calculation of the average recovery curve and its standard deviation be avoiding the memory-hungry vectorized approach and go for a iterative approach (two phased iterative approach will be necessary, i.e. one loop for the average curve and a second subsequent one for calculating the standard deviation). The real problem has something to do with the dask framework which shows memory leakage. Use memory profilinh tools to investigate this before changing RCC implementation (https://pythonhosted.org/Pympler/muppy.html)\n",
    "3. Bus error can be avoided by copying all auxilliary datasets (ranking database and motif annotation table) the node-local scratch storage ($VSC_SCRATCH_NODE).\n",
    "4. Incredibly slow even when using multiple cores. Do profiling to get more information on where the bottlenecks reside. 1 percent takes 30min so theoretically 60,000 modules would take 50 hours or 2 days.\n",
    "\n",
    "_Caveat:_ Still need to code conversion from regions back to genes to support AUcell in the last step of the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscenic.regions import RegionRankingDatabase, convert, Delineation, load\n",
    "from pyscenic.genesig import GeneSignature\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "DB_FOLDER = '/Users/bramvandesande/Projects/lcb/databases/'\n",
    "HG19_DB = 'hg19-regions-1M-9species.inverted.feather'\n",
    "RESOURCES_FOLDER = '/Users/bramvandesande/Projects/lcb/pyscenic/src/resources/tests'\n",
    "GMT_FNAME = 'c6.all.v6.1.symbols.gmt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dm6-5kb-upstream-full-tx-11species.mc8nr.feather\n",
      "dm6_UPSTREAM5KB_FULL_TX_motifRanking.RData\n",
      "hg19-500bp-upstream-10species.mc9nr.feather\n",
      "hg19-500bp-upstream-7species.mc9nr.feather\n",
      "hg19-regions-1M-9species.bed.gz\n",
      "hg19-regions-1M-9species.inverted.feather\n",
      "hg19-regions-1M-9species.inverted.identifiers.txt\n",
      "hg19-regions-220330-9species.extracted.feather\n",
      "hg19-regions-220330-9species.inverted.feather\n",
      "hg19-regions-220330-9species.inverted.identifiers.txt\n",
      "hg19-tss-centered-10kb-10species.mc9nr.feather\n",
      "hg19-tss-centered-10kb-7species.mc9nr.feather\n",
      "hg19-tss-centered-5kb-10species.mc9nr.feather\n",
      "hg19-tss-centered-5kb-7species.mc9nr.feather\n",
      "mm9-500bp-upstream-10species.mc9nr.feather\n",
      "mm9-500bp-upstream-7species.mc9nr.feather\n",
      "mm9-tss-centered-10kb-10species.mc9nr.feather\n",
      "mm9-tss-centered-10kb-7species.mc9nr.feather\n",
      "mm9-tss-centered-5kb-10species.mc9nr.feather\n",
      "mm9-tss-centered-5kb-7species.mc9nr.feather\n"
     ]
    }
   ],
   "source": [
    "!ls {DB_FOLDER}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "db = RegionRankingDatabase(os.path.join(DB_FOLDER, HG19_DB), name='hg19-1M-regions')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1223024"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(db.regions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "189"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "signatures = GeneSignature.from_gmt(os.path.join(RESOURCES_FOLDER, GMT_FNAME), 'HGNC', field_separator='\\t', gene_separator='\\t')\n",
    "len(signatures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.28067221,  8.68231655,  9.37239211,  9.22263887,  9.07035501,\n",
       "       13.21788222,  5.22955063, 11.14638916,  6.43318771,  3.91015555,\n",
       "       12.4239484 , 12.52129836, 13.09815328, 12.29158812,  3.6175327 ,\n",
       "       10.2431006 ,  3.25987139,  7.79902101,  6.3108416 ,  1.82386732,\n",
       "       13.28174089,  8.55206661,  2.6222048 ,  8.36171365, 11.23567103,\n",
       "        4.04145687, 10.40181721, 13.23567878,  1.96096473,  2.00766474,\n",
       "       10.2775063 ,  3.82094041, 11.08229546,  2.54086962,  8.04782619,\n",
       "        3.41068053,  8.87874084,  8.13214622,  2.79587443,  6.50719285,\n",
       "        5.19540032,  3.70436932,  9.47801332,  4.63939797,  8.91565459,\n",
       "        6.95054319,  7.07546034,  7.16014038, 11.29293814,  7.10672253])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.random.uniform(low=0.5, high=13.3, size=(50,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = signatures[89].copy(gene2weight=list(zip(signatures[89].genes, np.random.uniform(low=0.1, high=10.0, size=(len(signatures[89],))))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GeneSignature(name='PTEN_DN.V2_UP', nomenclature='HGNC', gene2weight=<frozendict {'IL32': 6.074815374172376, 'EHF': 9.001424216811124, 'BBC3': 4.2357212149387475, 'DENND2A': 6.054371626698191, 'CPA4': 6.3641123070290515, 'EMP1': 9.401430329587095, 'OSGIN1': 5.463787390838832, 'AKR1C3': 8.720942986125003, 'MYL10': 3.4105283810736533, 'IL17A': 4.812175424328964, 'GMIP': 3.6450138415728506, 'TRIB2': 6.96068973924421, 'BUD31': 3.582587321759326, 'LOC647835': 6.613040219776352, 'GPR87': 7.955599777815841, 'CCL16': 9.689740699757767, 'KRT15': 3.9560173119664617, 'NAV3': 7.925310025899901, 'LIPG': 9.768395717078745, 'ARVCF': 6.255898172533963, 'ZNF668': 0.3563703995239139, 'CEP41': 8.34206017733226, 'CLPB': 3.3049721274341075, 'HSPB8': 8.444547376459543, 'SERPINB5': 9.029959831725854, 'TSHB': 1.6675709830719216, 'ARHGDIB': 6.735016002274071, 'SATB1': 3.2187211360302883, 'ATP2B2': 2.7190791616140206, 'TMEM127': 8.436954307504214, 'SPOCK1': 7.244063009611832, 'LCN2': 3.779366793217826, 'ZNF549': 1.0897299324361671, 'ALDH1A3': 5.132920136403146, 'ERN2': 3.980506784717213, 'AGPAT4': 7.549262931407029, 'PRKY': 0.7669998996083734, 'PC': 6.157046477648348, 'ANGPTL4': 3.541297702323337, 'CUEDC1': 7.776393113134398, 'LAMA3': 0.3805861080250259, 'SAMD4A': 6.275959561266017, 'GNA14': 9.384582360460556, 'C14orf56': 6.592116895048843, 'NNAT': 3.000917248091861, 'RHBDF1': 3.832667356504324, 'CTNNBIP1': 1.4104002314335937, 'KCNJ3': 1.9535506091544175, 'KLHL28': 9.311591239103242, 'UPK2': 0.2553836687392682, 'COL13A1': 3.657801318334019, 'NID2': 1.2587778443354405, 'PCSK7': 1.3247124779049704, 'OTUB2': 6.415424957890448, 'KLK10': 1.9384483235204313, 'SERPINA5': 7.506983512515018, 'DEFB126': 7.096832377037246, 'BIRC3': 0.9548703523811214, 'SH3BGR': 5.935566841404641, 'NOL4': 4.006424895189592, 'CASP4': 9.5219517279095, 'LOXL2': 8.010245873607936, 'SLCO1B3': 9.301951786380243, 'IL4R': 9.748719111284682, 'PTGDS': 1.655710399370494, 'SPAM1': 5.2239274104138085, 'TP53I3': 1.6105944570030146, 'ANO2': 4.897798978612514, 'ZMAT3': 7.322107990047279, 'CLIP4': 7.868219960384064, 'ANXA3': 7.853786071185253, 'THBS4': 4.515415039885228, 'G0S2': 8.612165372724437, 'CHRDL1': 2.695738782506703, 'MIR22HG': 5.366801680650416, 'FYB': 4.642414907206539, 'PLAU': 7.299289272260676, 'OR12D2': 0.7200997515189347, 'GPM6A': 8.879236153341468, 'SPIB': 1.2370704612730044, 'CA5B': 3.0397126356585766, 'GK3P': 4.766886146791301, 'CHAC1': 8.30810256478099, 'C1orf20': 6.979875946177326, 'GVINP1': 4.31708937302265, 'SLITRK3': 8.815546350360673, 'KRT75': 4.891222825028978, 'IL1B': 8.54394468949546, 'CHRM4': 6.76635293943379, 'GLRA2': 5.170114792757705, 'TNFRSF10D': 4.072525866002369, 'JUN': 2.559685194552553, 'DGCR11': 1.5307678251918315, 'LIG4': 4.2905821960984465, 'GYS2': 6.318006944182711, 'DNAJB5': 8.3010361133243, 'UTRN': 9.123583468267775, 'OVGP1': 1.2463046802882733, 'GRID2': 7.4462129824513745, 'GPR183': 2.427451604591745, 'ST3GAL4': 8.367756076239164, 'PQLC3': 6.866673841764723, 'RND1': 0.38075139110313805, 'WT1-AS': 1.0035501398705997, 'PTGS2': 6.937643926478077, 'VAMP5': 0.6133480987891251, 'TBX5': 1.9867619959701737, 'VPS13C': 5.620782687619168, 'SLC2A4': 7.607554155799107, 'AHSP': 0.17435227167464987, 'NAALAD2': 9.267697568264095, 'SECTM1': 0.15475516826120617, 'CD86': 9.774313591611397, 'MPV17': 9.318126010743145, 'CEACAM21': 0.575660267740927, 'ALDH1A2': 9.687689520369126, 'IL8': 6.19223541171943, 'MON1B': 9.07305160066563, 'IL11RA': 1.7591152348593886, 'TNFRSF10C': 2.666649326987413, 'PPAP2B': 1.1476554265944177, 'ARL4C': 9.768582416278981, 'SLC1A3': 9.440107420390579, 'HSPB7': 4.51741148622951, 'ULK4': 4.755927557516935, 'APOL3': 8.099828136497669, 'HSPA4L': 2.1621195475754442, 'MC3R': 1.2864527826436472, 'KRT8P12': 6.097400324289782, 'ARHGEF5': 5.706118284840154, 'ANKMY1': 5.983154809796272, 'RTN1': 2.348677765237992, 'HTR1B': 2.6635572656882616, 'C12orf5': 6.526849405066924, 'BNC2': 3.7738179748497784, 'IL7R': 3.8067481039043516, 'AKR1B10': 9.818627300720786, 'ZMAT5': 2.940756077459487, 'KRT32': 9.253436799768224, 'P2RY2': 2.0790201094832974, 'LSP1': 0.9932325717197249, 'CSF1': 8.613291421560811, 'LOC440895': 4.0514221162669655}>)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GeneSignature(name='PTEN_DN.V2_UP', nomenclature='regions', gene2weight=<frozendict {'chr16-reg3999': 6.074815374172376, 'chr11-reg25202': 9.001424216811124, 'chr19-reg41736': 4.2357212149387475, 'chr19-reg41738': 4.2357212149387475, 'chr7-reg94938': 6.054371626698191, 'chr7-reg86082': 6.3641123070290515, 'chr12-reg12238': 9.401430329587095, 'chr16-reg59456': 5.463787390838832, 'chr16-reg59461': 5.463787390838832, 'chr10-reg4846': 8.720942986125003, 'chr7-reg68086': 3.4105283810736533, 'chr6-reg45817': 4.812175424328964, 'chr19-reg21924': 3.6450138415728506, 'chr2-reg11453': 6.96068973924421, 'chr7-reg65635': 3.582587321759326, 'chr3-reg103935': 7.955599777815841, 'chr17-reg29452': 9.689740699757767, 'chr17-reg34680': 3.9560173119664617, 'chr12-reg53504': 7.925310025899901, 'chr18-reg30389': 9.768395717078745, 'chr22-reg3047': 6.255898172533963, 'chr16-reg27291': 0.3563703995239139, 'chr16-reg27298': 0.3563703995239139, 'chr11-reg51268': 3.3049721274341075, 'chr12-reg84041': 8.444547376459543, 'chr18-reg40727': 9.029959831725854, 'chr1-reg97071': 1.6675709830719216, 'chr12-reg13517': 6.735016002274071, 'chr3-reg14907': 3.2187211360302883, 'chr3-reg14921': 3.2187211360302883, 'chr3-reg6822': 2.7190791616140206, 'chr3-reg6823': 2.7190791616140206, 'chr2-reg67224': 8.436954307504214, 'chr5-reg80584': 7.244063009611832, 'chr5-reg80585': 7.244063009611832, 'chr9-reg72570': 3.779366793217826, 'chr19-reg51859': 1.0897299324361671, 'chr15-reg64193': 5.132920136403146, 'chr16-reg21125': 3.980506784717213, 'chr6-reg109828': 7.549262931407029, 'chrY-reg681': 0.7669998996083734, 'chr11-reg45044': 6.157046477648348, 'chr11-reg45094': 6.157046477648348, 'chr19-reg9940': 3.541297702323337, 'chr17-reg49440': 7.776393113134398, 'chr18-reg13718': 0.3805861080250259, 'chr18-reg13938': 0.3805861080250259, 'chr14-reg19617': 6.275959561266017, 'chr14-reg19841': 6.275959561266017, 'chr9-reg30997': 9.384582360460556, 'chr20-reg27838': 3.000917248091861, 'chr16-reg78': 3.832667356504324, 'chr1-reg10626': 1.4104002314335937, 'chr2-reg105250': 1.9535506091544175, 'chr14-reg14166': 9.311591239103242, 'chr11-reg83107': 0.2553836687392682, 'chr10-reg45375': 3.657801318334019, 'chr14-reg17854': 1.2587778443354405, 'chr14-reg50521': 6.415424957890448, 'chr19-reg46157': 1.9384483235204313, 'chr14-reg51079': 7.506983512515018, 'chr20-reg19': 7.096832377037246, 'chr11-reg70051': 0.9548703523811214, 'chr21-reg17014': 5.935566841404641, 'chr21-reg17020': 5.935566841404641, 'chr18-reg19601': 4.006424895189592, 'chr18-reg19654': 4.006424895189592, 'chr18-reg19655': 4.006424895189592, 'chr11-reg71812': 9.5219517279095, 'chr11-reg71813': 9.5219517279095, 'chr11-reg71822': 9.5219517279095, 'chr8-reg16535': 8.010245873607936, 'chr12-reg16969': 9.301951786380243, 'chr16-reg23620': 9.748719111284682, 'chr9-reg83635': 1.655710399370494, 'chr7-reg81899': 5.2239274104138085, 'chr2-reg19508': 1.6105944570030146, 'chr12-reg6127': 4.897798978612514, 'chr3-reg120688': 7.322107990047279, 'chr2-reg24874': 7.868219960384064, 'chr4-reg45905': 7.853786071185253, 'chr5-reg47323': 4.515415039885228, 'chr1-reg152449': 8.612165372724437, 'chrX-reg40530': 2.695738782506703, 'chr5-reg24212': 4.642414907206539, 'chr10-reg49785': 7.299289272260676, 'chr6-reg24885': 0.7200997515189347, 'chr4-reg97690': 8.879236153341468, 'chr4-reg97760': 8.879236153341468, 'chr19-reg45463': 1.2370704612730044, 'chrX-reg6497': 3.0397126356585766, 'chr4-reg92422': 4.766886146791301, 'chr15-reg13136': 8.30810256478099, 'chr11-reg5971': 4.31708937302265, 'chr3-reg111900': 8.815546350360673, 'chr12-reg36108': 4.891222825028978, 'chr2-reg80146': 8.54394468949546, 'chr11-reg32543': 6.76635293943379, 'chrX-reg5864': 5.170114792757705, 'chr8-reg16253': 4.072525866002369, 'chr1-reg60087': 2.559685194552553, 'chr22-reg1859': 1.5307678251918315, 'chr13-reg52252': 4.2905821960984465, 'chr13-reg52256': 4.2905821960984465, 'chr12-reg17297': 6.318006944182711, 'chr9-reg19738': 8.3010361133243, 'chr9-reg19739': 8.3010361133243, 'chr6-reg95864': 9.123583468267775, 'chr1-reg93870': 1.2463046802882733, 'chr4-reg54026': 7.4462129824513745, 'chr13-reg46816': 2.427451604591745, 'chr11-reg90298': 8.367756076239164, 'chr2-reg9819': 6.866673841764723, 'chr12-reg32353': 0.38075139110313805, 'chr11-reg23236': 1.0035501398705997, 'chr1-reg135571': 6.937643926478077, 'chr2-reg63201': 0.6133480987891251, 'chr12-reg79917': 1.9867619959701737, 'chr15-reg28797': 5.620782687619168, 'chr17-reg7874': 7.607554155799107, 'chr16-reg27868': 0.17435227167464987, 'chr11-reg63829': 9.267697568264095, 'chr17-reg75055': 0.15475516826120617, 'chr3-reg80821': 9.774313591611397, 'chr3-reg80841': 9.774313591611397, 'chr2-reg22956': 9.318126010743145, 'chr19-reg36189': 0.575660267740927, 'chr15-reg25480': 9.687689520369126, 'chr15-reg25521': 9.687689520369126, 'chr4-reg42510': 6.19223541171943, 'chr16-reg53520': 9.07305160066563, 'chr9-reg19460': 1.7591152348593886, 'chr9-reg19462': 1.7591152348593886, 'chr8-reg16175': 2.666649326987413, 'chr1-reg58392': 1.1476554265944177, 'chr1-reg58393': 1.1476554265944177, 'chr2-reg161406': 9.768582416278981, 'chr5-reg22166': 9.440107420390579, 'chr5-reg22167': 9.440107420390579, 'chr5-reg22169': 9.440107420390579, 'chr1-reg16770': 4.51741148622951, 'chr3-reg29674': 4.755927557516935, 'chr22-reg18837': 8.099828136497669, 'chr22-reg18843': 8.099828136497669, 'chr22-reg18844': 8.099828136497669, 'chr4-reg72307': 2.1621195475754442, 'chr20-reg46288': 1.2864527826436472, 'chr7-reg97470': 5.706118284840154, 'chr2-reg168522': 5.983154809796272, 'chr14-reg23418': 2.348677765237992, 'chr14-reg23471': 2.348677765237992, 'chr14-reg23550': 2.348677765237992, 'chr6-reg56889': 2.6635572656882616, 'chr12-reg4625': 6.526849405066924, 'chr9-reg10133': 3.7738179748497784, 'chr5-reg21619': 3.8067481039043516, 'chr7-reg89932': 9.818627300720786, 'chr22-reg13060': 2.940756077459487, 'chr17-reg34601': 9.253436799768224, 'chr11-reg52145': 2.0790201094832974, 'chr11-reg2123': 0.9932325717197249, 'chr11-reg2136': 0.9932325717197249, 'chr1-reg92447': 8.613291421560811}>)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(test, db, Delineation.HG19_500BP_UP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5min 35s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r1 -n1\n",
    "db.load(convert(signatures[89], db, Delineation.HG19_500BP_UP))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6.4 (pyscenic_dev)",
   "language": "python",
   "name": "pyscenic_dev"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
